/*
 * Copyright (c) 2009, Luke Emrose
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * Redistributions of source code must retain the above copyright
 * notice, this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the
 * distribution.  Neither the name of Luke Emrose Sole Trader nor the
 * names of its contributors may be used to endorse or promote
 * products derived from this software without specific prior written
 * permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
 #ifndef __IMPROVED_PERLIN_H__
#define __IMPROVED_PERLIN_H__

#include "noisesGlobals.hpp"

/**
 * @file
 * Reference: http://mrl.nyu.edu/~perlin/noise/ - Java Reference Implementation Of Improved Noise - Copyright 2002 Ken Perlin
 * improved perlin noise
 *
 * Modified Code and sse, sse2 optimizations:
 * Luke Emrose - evolutionarytheory@gmail.com
 */

/**
 * @brief
 * permutation table
 */
static int permutation[] =
{
	151, 160, 137, 91, 90, 15, 
	131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 
	190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 
	88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168,  68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 
	77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 
	102, 143, 54,  65, 25, 63, 161,  1, 216, 80, 73, 209, 76, 132, 187, 208,  89, 18, 169, 200, 196, 
	135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186,  3, 64, 52, 217, 226, 250, 124, 123, 
	5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 
	223, 183, 170, 213, 119, 248, 152,  2, 44, 154, 163,  70, 221, 153, 101, 155, 167,  43, 172, 9, 
	129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 
	251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 
	49, 192, 214,  31, 181, 199, 106, 157, 184,  84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 
	138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
};

/**
 * @brief
 * We use a per-frame initalizer to create the table when the plugin is 
 * first loaded in the renderer. Since this table is read only, we do not
 * need to worry about locking multi-thread access to this global variable.
 */
static int p[ 512 ];

/**
 * @brief
 * precalculate noise table
 */
void initImprovedPerlinNoiseTable( RixContext *ctx )
{
	for( int i = 0; i < 256; ++i )
	{
		p[ 256 + i ] = p[ i ] = permutation[ i ];
	}
}

/**
 * @brief
 * smooth step function which has zero first and second derivatives at both t = 0 and 1
 */
template < typename T >
inline static T fade( const T t ) 
{
	return ( t * t * t * ( t * ( t * 6 - 15 ) + 10 ) );
}

#ifdef USE_SSE
/**
 * @brief
 * sse optimized smooth step function which has zero first and second derivatives at both t = 0 and 1
 */
template <>
inline static __m128 fade( const __m128 t ) 
{
	// perform t * 6
	const __m128 r0 = _mm_mul_ps( t, _mm_set1_ps( 6.0f ) );

	// perform [ t * 6 ] - 15
	const __m128 r1 = _mm_sub_ps( r0, _mm_set1_ps( 15.0f ) );

	// perform t * [ t * 6 - 15 ]
	const __m128 r2 = _mm_mul_ps( t, r1 );

	// perform [ t * 6 - 15 ] + 10
	const __m128 r3 = _mm_add_ps( r2, _mm_set1_ps( 10.0f ) );

	// perform t * t * t * [ ( t * 6 - 15 ) + 10 ]
	__m128 r4 = _mm_mul_ps( r3, t );
	r4 = _mm_mul_ps( r4, t );
	r4 = _mm_mul_ps( r4, t );

	// obtain result
	return r4;
}
#endif

/**
 * @brief
 * linear interpolation
 */
template < typename T >
inline static T lerp( const T t, const T a, const T b ) 
{
	return ( a + t * ( b - a ) ); 
}

#ifdef USE_SSE
/**
 * @brief
 * sse optimized linear interpolation
 */
template <>
inline static __m128 lerp( const __m128 t, const __m128 a, const __m128 b ) 
{
	// b - a
	const __m128 r0 = _mm_sub_ps( b, a );

	// t * [ ( b - a ) ]
	const __m128 r1 = _mm_mul_ps( t, r0 );

	// a + [ t * ( b - a ) ]
	const __m128 r2 = _mm_add_ps( a, r1 );

	// obtain result
	return r2;
}
#endif

/**
 * @brief
 * gradient calculation
 */
template < typename T >
inline static T grad( const int hash, const T x, const T y, const T z ) 
{
	int h = hash & 15;                   // CONVERT LO 4 BITS OF HASH CODE
	T u = h < 8 ? x : y;                 // INTO 12 GRADIENT DIRECTIONS.
	T v = h < 4 ? y : h == 12 || h == 14 ? x : z;
	return ( ( ( h & 1 ) == 0 ? u : -u ) + ( ( h & 2 ) == 0 ? v : -v ) );
}

/**
 * @brief
 * noise calculation
 */
template < typename T, typename U >
inline static T improvedNoise_impl( const U in_pt ) 
{
	int X, Y, Z;
	int A, B, AA, AB, BA, BB;

	// calculate floor of values
	U f;
	f[ 0 ] = floor( in_pt[ 0 ] );
	f[ 1 ] = floor( in_pt[ 1 ] );
	f[ 2 ] = floor( in_pt[ 2 ] );

	// find unit cube that contains point
	X = ( ( int ) f[ 0 ] ) & 255;	
	Y = ( ( int ) f[ 1 ] ) & 255;	
	Z = ( ( int ) f[ 2 ] ) & 255;

	// find relative x, y, z of point in cube
	U pp;
	pp[ 0 ] = in_pt[ 0 ] - f[ 0 ];
	pp[ 1 ] = in_pt[ 1 ] - f[ 1 ];
	pp[ 2 ] = in_pt[ 2 ] - f[ 2 ];

	// compute fade curves for each of x, y, z
	U uvw;
	uvw[ 0 ] = fade( pp[ 0 ] );
	uvw[ 1 ] = fade( pp[ 1 ] );
	uvw[ 2 ] = fade( pp[ 2 ] );

	// hash coordinates of the 8 cube corners
	A = p[ X ] + Y;
	AA = p[ A ] + Z;
	AB = p[ A + 1 ] + Z;
	B = p[ X + 1 ] + Y;
	BA = p[ B ] + Z;
	BB = p[ B + 1 ] + Z;

	// and add blended results from 8 corners of cube
	return	lerp( uvw[ 2 ],	lerp( uvw[ 1 ],	lerp( uvw[ 0 ],	grad( p[ AA ],		pp[ 0 ],		pp[ 1 ],		pp[ 2 ] ),
															grad( p[ BA ],		pp[ 0 ] - 1,	pp[ 1 ],		pp[ 2 ] ) ),
											lerp( uvw[ 0 ],	grad( p[ AB ],		pp[ 0 ],		pp[ 1 ] - 1,	pp[ 2 ] ),
															grad( p[ BB ],		pp[ 0 ] - 1,	pp[ 1 ] - 1,	pp[ 2 ] ) ) ),
							lerp( uvw[ 1 ],	lerp( uvw[ 0 ],	grad( p[ AA + 1 ],	pp[ 0 ],		pp[ 1 ],		pp[ 2 ] - 1 ),
															grad( p[ BA + 1 ],	pp[ 0 ] - 1,	pp[ 1 ],		pp[ 2 ] - 1 ) ),
											lerp( uvw[ 0 ],	grad( p[ AB + 1],	pp[ 0 ],		pp[ 1 ] - 1,	pp[ 2 ] - 1 ),
															grad( p[ BB + 1 ],	pp[ 0 ] - 1,	pp[ 1 ] - 1,	pp[ 2 ] - 1 ) ) ) );
}

#ifdef USE_SSE
/**
 * @brief
 * sse optimized noise calculation
 */
template < typename T >
inline static T improvedNoise_impl( const __m128 in_pt ) 
{
	// calculate floor of values
	const __m128 r0 = _mm_floor_ps( in_pt );

	// find relative x, y, z of point in cube
	const __m128 r1 = _mm_sub_ps( in_pt, r0 );

	// first convert values to integers
	const __m128i r2 = _mm_cvttps_epi32( r0 );

	// find unit cube that contains point
	const __m128i r3 = _mm_and_si128( r2, _mm_set1_epi32( 255 ) );

	// compute fade curves for each of x, y, z
	const __m128 r4 = fade( r1 );

	// hash coordinates of the 8 cube corners

	// extract X, Y, Z
	int X = ( ( int* ) &r3 )[ 0 ];
	int Y = ( ( int* ) &r3 )[ 1 ];
	int Z = ( ( int* ) &r3 )[ 2 ];

	// calculate A's, B's
	int A = p[ X ] + Y;
	int AA = p[ A ] + Z;
	int AB = p[ A + 1 ] + Z;
	int B = p[ X + 1 ] + Y;
	int BA = p[ B ] + Z;
	int BB = p[ B + 1 ] + Z;

	// precalculations of lookups
	int pAA = p[ AA ];
	int pBA = p[ BA ];
	int pAB = p[ AB ];
	int pBB = p[ BB ];

	// precalculations of lookups for + 1
	int pAA_1 = p[ AA + 1 ];
	int pBA_1 = p[ BA + 1 ];
	int pAB_1 = p[ AB + 1 ];
	int pBB_1 = p[ BB + 1 ];

	// get x, y, z
	_CRT_ALIGN( 16 ) float xyz[ 4 ];
	_mm_store_ps( xyz, r1 );
	float x = xyz[ 0 ];
	float y = xyz[ 1 ];
	float z = xyz[ 2 ];

	// get u, v, w
	_CRT_ALIGN( 16 ) float uvw[ 4 ];
	_mm_store_ps( uvw, r4 );
	float u = uvw[ 0 ];
	float v = uvw[ 1 ];
	float w = uvw[ 2 ];

	// perform all grads
	float gAA = grad( pAA, x    , y    , z );
	float gBA = grad( pBA, x - 1, y    , z );
	float gAB = grad( pAB, x    , y - 1, z );
	float gBB = grad( pBB, x - 1, y - 1, z );

	float gAA_1 = grad( pAA_1, x    , y    , z - 1 );
	float gBA_1 = grad( pBA_1, x - 1, y    , z - 1 );
	float gAB_1 = grad( pAB_1, x    , y - 1, z - 1 );
	float gBB_1 = grad( pBB_1, x - 1, y - 1, z - 1 );

	// perform lerps
	float u_gAA_gBA = lerp( u, gAA, gBA );
	float u_gAB_gBB = lerp( u, gAB, gBB );

	float u_gAA_1_gBA_1 = lerp( u, gAA_1, gBA_1 );
	float u_gAB_1_gBB_1 = lerp( u, gAB_1, gBB_1 );

	float v_u_gAA_gBA_u_gAB_gBB = lerp( v, u_gAA_gBA, u_gAB_gBB );
	float v_u_gAA_1_gBA_1_u_gAB_1_gBB_1 = lerp( v, u_gAA_1_gBA_1, u_gAB_1_gBB_1 );

	// and add blended results from 8 corners of cube
	return lerp( w, v_u_gAA_gBA_u_gAB_gBB, v_u_gAA_1_gBA_1_u_gAB_1_gBB_1 );
}
#endif

template < typename T, typename U >
inline static T improvedNoiseFBM_impl( const U in_pt, const T numOctaves, const T lacunarity, const T gain )
{
	// initialize all octave variables
	T amp = 1;
	T freq = 1;
	T norm = 0;
	T noiseSum = 0;

	U pt;
	pt[ 0 ] = in_pt[ 0 ];
	pt[ 1 ] = in_pt[ 1 ];
	pt[ 2 ] = in_pt[ 2 ];

	// iterate over all octaves, and sum up each octave of noise
	for( int octave = 0; octave < numOctaves; ++octave )
	{
		// calculate noise
		pt[ 0 ] *= freq;
		pt[ 1 ] *= freq;
		pt[ 2 ] *= freq;
		noiseSum += amp * improvedNoise_impl< T >( pt );

		// update values for next iteration
		norm += amp;
		freq *= lacunarity;
		amp *= gain;
	}

	return noiseSum / norm;
}

#ifdef USE_SSE
template < typename T >
inline static float improvedNoiseFBM_impl( const __m128 in_pt, const T numOctaves, const T lacunarity, const T gain )
{
	// copy in the input point for manipulation, to ensure that this function has no side effects
	__m128 pt = in_pt;

	// initialize all octave variables ( noiseSum, norm, freq, amp )
	// remember that these go into the register in reverse order!
	__m128 vals = _mm_set_ps( 1, 1, 0, 0 );

	// set the constant multiplier for gain and lacunarity
	const __m128 gain_lacunarity = _mm_set_ps( gain, lacunarity, 1, 1 );

	// iterate over all octaves, and sum up each octave of noise
	for( int octave = 0; octave < numOctaves; ++octave )
	{
		// calculate noise
		// perform pt *= freq
		pt = _mm_mul_ps( pt, _mm_set1_ps( ( ( float* ) &vals )[ 2 ] ) );

		// perform:
		// noiseSum += amp * improvedNoise_impl< T >( pt );
		// and
		// norm += amp;
		// simultaneously
		vals = _mm_add_ps( vals, _mm_set_ss( ( ( float* ) &vals )[ 3 ] * improvedNoise_impl< T >( pt ) ) );
		vals = _mm_add_ps( vals, _mm_set_ps( 0, 0, ( ( float* ) &vals )[ 3 ], 0 ) );

		// perform:
		// freq *= lacunarity
		// and
		// amp *= gain
		// simultaneously
		vals = _mm_mul_ps( vals, gain_lacunarity );
	}

	//return noiseSum / norm;
	return ( ( float* ) &vals )[ 0 ] / ( ( float* ) &vals )[ 1 ];
}
#endif

template < typename T, typename U >
inline static T improvedNoiseTurbulence_impl( const U in_pt, const T numOctaves, const T lacunarity, const T gain )
{
	// initialize all octave variables
	T amp = 1;
	T freq = 1;
	T norm = 0;
	T noiseSum = 0;

	U pt;
	pt[ 0 ] = in_pt[ 0 ];
	pt[ 1 ] = in_pt[ 1 ];
	pt[ 2 ] = in_pt[ 2 ];

	// iterate over all octaves, and sum up each octave of noise
	for( int octave = 0; octave < numOctaves; ++octave )
	{
		// calculate noise
		pt[ 0 ] *= freq;
		pt[ 1 ] *= freq;
		pt[ 2 ] *= freq;
		noiseSum += amp * fabs( improvedNoise_impl< T >( pt ) );

		// update values for next iteration
		norm += amp;
		freq *= lacunarity;
		amp *= gain;
	}

	return noiseSum / norm;
}

#endif